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The relative sidereal variation in the arrival direction of primary cosmic ray nuclei of median energy 
10 TeV was measured using downward, through-going muons detected with the Super-Kamiokande-I 
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detector. The projection of the anisotropy map onto the right ascension axis has a first harmonic 
amplitude of (6.64 ±0.98 stat. ±0.55 syst.) x 10~* and a phase at maximum at (33.2° ±8.2° stat. ± 
5.1° syst.) right ascension. A sky map in equatorial coordinates indicates an excess region in 
the constellation of Taurus and a deficit region toward Virgo. The excess region is centered at 
{ot, St) = (75° ± 7°, -5° ± 9°) with a half opening angle xt = (39 ± 7)°; the excess flux is 
(0.104 ± 0.020)% above the isotropic expectation. The corresponding parameters for the deficit 
region are (ay, Sy) = (205° ± 7°, 5° ± 10°), Xv = (54 ± 7)°, and (-0.094 ± 0.014)%. The data do 
not allow us to rule out a pure dipole form for the anisotropy (allowed at 13% confidence level); they 
are better described by the excess and deficit cones described above. We explored the implications 
under the assumption that the true anisotropy is not distorted too much by the analysis filter so 
that it is well-described by the observed excess and deficit cones. 

PACS numbers: 95.85.Ry, 96.50.Bh 



I. INTRODUCTION 

The flux of cosmic rays with energy per nucleon in 
the range 10^^ ~ 10^"* eV is known to have a sidereal 
anisotropy of several times 10^'*. The anisotropy is due 
to a combination of effects. Conipton and Getting ^ pro- 
posed in 1935 that the motion of the solar system relative 
to the rest frame of the cosmic ray plasma should cause an 
energy-independent dipole anisotropy whose maximum is 
in the direction of this motion. Solar diurnal and sea- 
sonal changes in the atmospheric temperature can in- 
duce a sidereal variation in the cosmic ray rate 0]. The 
anisotropy that remains after accounting for these effects 
is presumably of Galactic origin, with possible modula- 
tions due to the hehosphere (see, for example, [1,|3]) and, 
at the lowest energies, solar wind and magnetic field. 

In this article, we present a report on the observation 
of cosmic ray anisotropy with the Super-Kamiokande I 
(SK-I) detector. SK-I can make a unique contribution to 
this subject because of the large overburden and detector 
size. The overburden makes SK-I sensitive to primary 
cosmic ray energies normally attainable with extensive 
air shower arrays, while the large statistics and excellent 
muon tracking resolution enabled the creation of a two 
dimensional map of the anisotropy, which is the first- 
published muon-bascd map Q. 



II. THE DETECTOR AND THE DATA 

SK-I is a 50 kiloton underground imaging water 
Cherenkov detector in Kamioka, Japan at geographical 
coordinates 36°25'32.6" N, 137°18'37.1" E and an alti- 
tude of 370 m above sea level. The vertical overburden 
is about 1000 meters, or 2700 meters water equivalent. 
The detector's design was optimized for the detection of 
neutrinos and nucleon decay; the placing of the detec- 
tor under large overburden to shield against cosmic ray 
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muons is an important part of this design. The over- 
burden shields all charged cosmic ray secondaries except 
muons with energy above 0.8 TeV. The portion of the de- 
tector sensitive to muons is a cylinder of diameter 33.8 m 
and height 36.2 m, giving a target area between 1000 m^ 
and 1200 m^ depending on the zenith angle. The aver- 
age cosmic ray muon event rate was 1.8 Hz. More details 
about the SK-I detector are reported in Q. 

The data used in this analysis were collected between 
June 1, 1996 and May 31, 2001. During this period, the 
detector was live for 1662.0 days (91.0% of the time) and 
registered 2.54 x 10® muon events. Muon track recon- 
struction was performed with an algorithm developed in 
SK-I to examine the spatial correlation between spalla- 
tion products and parent muons in the solar neutrino 
analysis [7|. In order to maintain the angular resolution 
within 2°, the muon tracks were required to be longer 
than 10 m and downward-going; 82.6% of the events 
(2.10 X 10®) satisfied these requirements. The reliabil- 
ity of using muon tracks for astronomical purposes was 
confirmed by the observation of the shadow of the moon 
and the sun Q. 

The relationship between the energy of the detected 
muon and the energy of the primary cosmic ray that 
produced it is described by a response function (see, for 
example, Ref. @). For SK-I, the threshold muon en- 
ergy is 0.8 TeV for the thinnest part of the overburden. 
The corresponding median primary cosmic ray energy is 
about ten times larger [Q|, while the spread in the pri- 
mary cosmic ray energy is about an order of magnitude 
above and below the median. SK-I is, therefore, sen- 
sitive to primary cosmic rays with energy in the range 
several TeV to several hundred TeV. 

The variation in the overburden along different lines 
of sight explains most of the features seen in the muon 
event rate in the horizontal coordinate system (Fig. [1]) . 
This variation implies that the muon threshold energy, 
and, therefore, the median primary cosmic ray energy, 
vary with direction. A given point in the celestial coor- 
dinate system traces out a trajectory of fixed declination 
as the Earth rotates; the overburden along the line of 
sight to this point varies with this motion. For instance, 
the thickness of the overburden at the apex of the decli- 
nation = 0° trajectory is about 2300 meters water equiv- 
alent, which corresponds to a median primary cosmic 
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FIG. 1: Event rate (day~^m~'^sr~^) in horizontal coordinates. 
The dotted curves indicate contours of constant declination, 
while the arrows indicate the apparent motion of stars with 
the rotation of the earth. 

ray energy of about 8 TeV. In contrast, the correspond- 
ing thickness and median energy for declination = 50° S 
are 3800 meters water equivalent and 20 TeV, respec- 
tively. The average overburden over one rotation period 
is, therefore, a function of declination, which implies that 
the primary cosmic ray spectrum seen by SK varies with 
declination. For this reason, the anisotropy presented 
here is a spectrum-weighted average over a broad range 
of energies spanning several to several hundred TeV. 

It is seen from Fig. [T] that the detector is most sensi- 
tive to cosmic rays originating from the south. During 
5 years' operation of SK-I, the most sensitive direction 
was exposed to every 1° RA slice of the celestial sphere 
for an average of about 4.6 days. The exposure, however, 
was unequal for different directions because the detector 
was dead for several minutes almost every day, with occa- 
sional periods of down time lasting many hours. Periods 
of detector down time were not random, but occurred 
more often during day time hours. These periods of dead 
time, accumulated over five years, introduce fluctuations 
in the exposure times for different directions. The maxi- 
mum and minimum deviations from the average exposure 
were ■ In the celestial coordinate analysis of Sec. lIIIi 
the rates were corrected to account for different measured 
detector live time in different sidereal time bins. 

The atmosphere is a part of the detector in the sense 
that it is responsible for converting primary cosmic 
rays into muons that can penetrate the overburden. It 
is a dynamic detector component because its density 
changes with temperature and pressure, and the muon 
rate changes accordingly. The relative variation in the 
muon rate due to atmospheric variation has relatively 
strong Fourier components with frequencies correspond- 
ing to one year and one solar day. The solar diurnal com- 
ponent of the muon rate is, to some extent, modulated by 
a seasonally varying signal, giving rise to spurious side- 
real variation. The magnitude of this variation is large 
(comparable to the true magnitude) when the observa- 
tion period is restricted seasonally. However, averaged 



over exact-year periods, the spurious variation largely 
cancels out; we estimate it to be only 18% of the true 
magnitude. For this reason, we chose our data to span 
the exact five year period between June 1, 1996 and June 
1, 2001. In the celestial coordinate analysis below, we 
statistically subtract this spurious variation from the ob- 
served signal using the method of Farley and Storey 0] . 
This subtraction introduces an uncertainty of about 10% 
to the true magnitude of the sidereal anisotropy [35j . 
Technical details on the subtraction of the spurious atmo- 
spheric effect is given in Appendix]^ It is also important 
to verify that the input to this correction is sound; we 
show in Appendix [B] that the monthly muon rate vari- 
ation, which is primarily due to atmospheric effects, is 
properly measured in Super-Kamiokande. 

Another correction that can be made to the anisotropy 
measurement is that for the Compton-Getting effect. We 
chose not to subtract this in our main result because the 
rest frame of the cosmic ray plasma - an important in- 
put for the subtraction - is not known. In principle this 
can be measured with our data by measuring the sea- 
sonal change in the anisotropy. When the earth's orbital 
phase is such that the orbital velocity moves against the 
bulk cosmic ray motion, the flux is enhanced in this di- 
rection; six months later, the effect should be smaller, 
as the earth moves along with the bulk motion and the 
relative velocity between the observer and the cosmic 
ray rest frame is at a minimum. In practice, however, 
the cosmic ray rest frame was not measurable in this 
manner because the seasonal change in the flux intro- 
duced by atmospheric effects was much larger than the 
Compton-Getting signal. In the literature, two cosmic 
ray rest frames are often assumed: (1) the local stan- 
dard of rest (the frame of the average motion of stars in 
the neighborhood of the sun); (2) the rest frame of the 
local interstellar medium. Appendix [C] gives the result 
of anisotropy measurements with the subtraction of the 
Compton-Getting effect assuming these two rest frames, 
and the resulting anisotropy parameters are summarized 
in tables [IV] and El 



III. CELESTIAL COORDINATE ANALYSIS 

In this section, we describe the celestial coordinate 
analysis of the anisotropy. First, the general mathemat- 
ical framework of the analysis is given. An important 
part of this discussion is the distortion introduced to the 
anisotropy map by the analysis method. In the second 
subsection, the implementation of the method to the data 
is described. In the third subsection, the anisotropy in 
one dimension (i.e. as a function of right ascension) is 
presented. This is followed by the anisotropy map in two 
dimensions. 
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A. Analysis Method: Mathematical Framework 
and the Distortion Introduced to the Anisotropy 
Map 

The number of downward-going muon events from the 
celestial coordinates (a, S) observed in SK-I may be ex- 
pressed as follows: 



E 

k 



[1 + e{ai,Sj)] ■ R{ai - Tk,5j) ■ Tk ■ dVti 



(1) 

The indexes i, j, and k are for the right ascension, decli- 
nation, and sidereal time, Nsid is the number of sidereal 
time bins, and is the total detector live time in the 
sidereal time bin k. The function R{a — r, 5) is the dif- 
ferential event rate for an isotropic cosmic ray flux in the 
detector coordinate system (parameterized by hour angle 
and declination) . The form of this function is determined 
by the overburden and the zenith angle dependence of the 
cosmic ray flux. The function e(a, 5) represents the true 
cosmic ray anisotropy. 

As stated in the previous section, the term Tj in Eqn.[T] 

varies by !'^}'3%- This variation of purely instrumental 
origin was removed by multiplying each sidereal time bin 
by the following weight: 



Wk 



Tk 



(2) 



Once this correction is made, the raw anisotropy data 
can be expressed as follows: 

n{a.„5j) = '^[l + e{ai,5j)]- R{ai-Tk,5j) (3) 

k 

= [\ + e{a,,5j)]-p{5j) 

In the second line, p{5j) — J2k"'' ^i^i ~ '^k,Sj) - i.e. 
the summation of hour angle erases the right ascension 
dependence because of the cyclical nature of the function 
R. In other words, the sum is independent of the starting 
point, specified by a^. 

Ideally, one would like to extract the anisotropy func- 
tion e(a, S) from the data. In practice, this cannot be 
done because the declination dependence for an isotropic 
flux, p{S), can neither be measured from the data nor cal- 
culated to an accuracy required to extract an anisotropy 
of order several parts per 10,000. However, p{5) can be 
factored out by calculating the following ratio: 



1 

[1 + 



i 



(6) 



(7) 



The second term in Eqn. [5] distorts the true anisotropy 
function e(a, i5) in an unknown, but restricted way. As an 
illustration of the nature of the distortion, let us imagine 
that, for a fixed declination S, e(a, S) is well-described 
by a sinusoidal function; this function can be written as 
a sum of a constant offset and a sinusoidal term 

whose average is zero. The second term in Eqn.[5]removes 
the constant offset. In more precise mathematical terms, 
this distortion can be described as follows. The spherical 
harmonic decomposition of e(a, 5) is: 



(8) 



The angle A = 7r/2 — (5 is the complement of the decli- 
nation, and it is measured relative to the z axis (Earth's 
rotation axis) in the usual notation for spherical harmon- 
ics. The modified anisotropy function A(a, 6) is related 
to this as follows: 



A{a,S) = aj^, 



r,,™(a,A)-— / da' Ye^ra{a',X) 
Zn 



(9) 



The new coefficients bg^m have the following values: 



m = 

ai,m m ^ 



(10) 



It is seen that the axisymmetric terms (i.e. m — terms) 
are zeroed out. 

As a concrete example, consider the effect of the distor- 
tion on the first harmonic of an axisymmetric anisotropy 
(i.e. a dipole anisotropy) of magnitude D along an ar- 
bitrary direction (q!o,(5o) in equatorial coordinates. The 
anisotropy function has the following form: 



A{ai,Sj 



eiaz,Sj) - {e{6j)) 



(4) 
(5) 



The second line is an approximation that ignores second 
and higher order terms in e. The averages indicated by 
the brackets is over right ascension bins. Explicitly, 



A{a, S) ^ D [cos 5 cos Sq cos(a — uq) + sin ^ sin Sq] (11) 
After subtracting the constant offset, this becomes: 

A{a,S) = I?cos(5ocos<5cos(a — ao) (12) 
= I)cos(5cos(a — ao) (13) 
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The second line is the form of the anisotropy for the 
projection of the original dipole in the equatorial plane; 
the dipole strength D = D cos Sq is the length of this 
projection. 

Functions with higher harmonics behave in more com- 
plicated ways. For instance, consider an anisotropy func- 
tion that can be described with two cones, one with 
excess flux, the other with a deficit flux. If these two 
cones are not 180° opposite one another, significant con- 
tributions from higher harmonics must, by necessity, be 
present. If the declination of the two cones is similar, 
then the distortion is small (the excess cone cancels out 
the deficit cone, making the constant offset small). If 
they are different, the distortion pushes the declination 
of the two cones towards each other, while the right as- 
cension is not affected at all. 



B. Analysis Method: Implementation 

In the analysis presented here, the cosmic ray 
anisotropy in the celestial sphere was plotted in one and 
two dimensions. The two dimensional map corresponds 
to A{a, 6) shown in Eqn. [4l The one dimensional map 
can be plotted either as a function of right ascension or 
sidereal time. In the former, the plot corresponds to the 
following: 



a{ai) 

where m{a) is defined as: 



with the correction e about 1%). The relative variation of 
this histogram about its mean corresponds to a(a). The 
two dimensional anisotropy A{a, 5) is made exactly like 
a(a), but in 10° strips of declination. Finally, the one 
dimensional plot a(T) is made by making a histogram of 
the number of muon events in sidereal time bins, dividing 
this bin-by-bin with a histogram of the detector live time 
in sidereal time bins, and taking the variation relative to 
the mean. The function a(r) can also be thought of as 
a(a) with a replaced with t. The shape of the resulting 
function is similar to that of a{a) because R{a — r, S) is 
generally largest when a ft! r (Fig.[T]). In other words, at 
any given moment, the right ascension of a muon track 
measured with SK is approximately equal to the side- 
real time, or to the right ascension of the zenith. For 
this reason, we shall refer to a^r) as the 'zenith-type' 
anisotropy, while a{a) shall be referred to as the 'track- 
type' anisotropy because it is made using information 
from muon tracks. The 'zenith-type' anisotropy is equiv- 
alent to smearing the 'track-type' anisotropy. Clearly, the 
function a{a) is a better probe of cosmic ray anisotropy 
than a(r), but we have, nevertheless, produced a(T) be- 
cause most other underground muon measurements are 
presented in this way. 

Spurious sidereal variation of atmospheric origin de- 
scribed in Sec.lTTlwas subtracted from all plots and maps 
unless otherwise noted. The spurious variation has little 
effect on the best fit value of the parameters describing 
the anisotropy, but it significantly increases the uncer- 
tainty. 



C. Right Ascension Distribution 



(15) 



(m) corresponds to m{a) averaged over right ascension, 
and the function n{a, 6) is defined in Eqn. [3) The one 
dimensional map as a function of sidereal time is defined 
as follows: 



A track-type plot of the right ascension of cosmic rays 
before subtracting the spurious sidereal anisotropy from 
atmospheric effects is shown as data points in Fig. [3 (a). 
The solid curve is the best fit of the first two harmonics 
to the data, while the dashed curve (almost overlapping 
the solid one) is the sidereal variation after correcting for 
the atmospheric effect. The curves are parameterized as 
follows: 



a(Tfc) 



m{Tk) - (m) 
(m) 



The function rh{Tk) is defined as: 



(16) 



No, Ns 

"rniTk) = y^y^ [1 + ^{ai,5j)] ■ R{ai - Tk,6j), (17) 



and (m) is the average of m(T) over sidereal time. 

In practice, the one dimensional anisotropy plot a{a) 
is made by first making a histogram of the muon track 
right ascension, where each entry is weighted by Wk in 
Eqn. [5] in order to equalize the exposure to all directions 
in the celestial sphere (the value of the weight is 1 ± e, 



F{x) = Ai-cos 



180 ^ 



[I8O ^ ^ ' 

(18) 

The best fit parameters are summarized in Table HI 
The parameter errors are statistical, except measure- 
ments with atmospheric correction (track/corr. and 
zenith/corr.), where the first error is statistical and 
the second error is the systematic error introduced by 
subtraction of the atmospheric effect (see Appendix [X)) . 

The zenith-type plot of cosmic ray right ascension is 
shown in Fig.[2](b). The fit parameters of (b) are similar 
to those of (a), but the amplitudes are smaller, which is 
consistent with the fact that (b) is obtained by smearing 
(a). 
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1.002 



1.001 




0.999 



180 240 300 360 
Right Ascension (deg) 



1.001 



0.999 - 



0.991 




180 240 300 360 

Right Ascension at Zenith (deg) 

FIG. 2: (a) Track-type right ascension projection plot, (b) 
Zenith- type plot. The error bars represent statistical errors. 
The solid curve in each frame is the best fit of the first two 
harmonic functions. The dashed curve (almost overlapping 
the solid curve) is the first two harmonics after subtracting 
the atmospheric contribution. 



Anisotropy measurements based on zenith-type plots 
from Kamiokande [lOl and MACRO (both deep un- 
derground experiments like SK-I) are also summarized in 
Table [H Only the first ter m in Eqn. [18] was fit to their 
data. The amplitude and phase show good agreement 
with those of SK-I. 

Figure [3] shows the amplitude and phase of the best- fit 
first harmonic function fit to zenith-type plots from nu- 
merous experiments. The SK-I result is consistent with 
the trend. 



D. Sky Map of the Anisotropy 

An anisotropy map in the celestial sphere is obtained 
by making track-type plots in 10° strips of declination, 
giving 10° X 10° pixels. The result is shown in Fig. [4] 
(smoothing applied for visual purposes only). Each 
10° X 10° cell in map (a) shows the fractional variation 
from the isotropic case, while that in (b) shows the stan- 
dard deviation of this variation. The fractional varia- 
tions from isotropy are large and erratic in (a) because 
of oversampling, but the positive and negative variations 
are clearly clustered. In order to optimize the binning, as 
well as to identify and characterize the excess and deficit 
regions, we used a clustering algorithm described in Ap- 
pendix[D] The algorithm indicates a unique region of ex- 
cess in the direction of the constellation Taurus (ar, St) 



= (75°±7°, — 5°±9°), while a unique region of deficit was 
found in the constellation of Virgo {ay, Sy) = (205° ±7°, 
5° ± 10°). The half opening angle of the "Taurus" region 
is 39° ± 7° with a relative rate (0.104 ± 0.020)% above 
average, while the size of the "Virgo" region is 54° ± 7° 
with a relative rate of (0.094 ± 0.014)% below average. 
All errors are statistical (see Appendix [Pl for the method 
used to obtain the statistical error of reconstructed cone 
parameters). These results are summarized in Table HIl 

The observed anisotropy is unlikely to be due to a ran- 
dom fluctuation of an isotropic cosmic ray flux. The cal- 
culation of the statistical significance and the rejection of 
the null hypothesis are performed as follows. The number 
of events in the Taurus excess cone is (0.104 ± 0.020)% 
above the expectation from the isotropic distribution, 
which corresponds to a gaussian probability of 2.0 x 10^ 
However, since the entire sky above the horizon was 
searched with a variable half-opening angle, the actual 
probability for this sort of deviation to occur is larger by 
some trials factor. In order to determine this, 1 x 10^ 
isotropic sky maps were generated with statistical fluc- 
tuations generated with a random number generator. To 
cover the angular size of the Taurus excess, we counted 
the number of maps with reconstructed cone with in- 
cone standard deviation > 5.2 sigma and half-opening 
angle between 30° and 60°. The number of such maps 
was 378 out of 10^ generated maps, giving a post-trials 
probability of 3.78 x 10^^. Similarly, the number of 
events in the Virgo deficit cone is (0.094 ±0.014)% below 
the isotropic expectation, which is a 6.7 standard de- 
viation effect corresponding to a gaussian probability of 
2.1 X 10~^^. Among the 1 x 10^ generated maps, none had 
a deviation as large as observed with half-angle between 
30° and 60°. We, therefore, set a 90% confidence level 
upper limit of the post-trials probability at 2.3 x 10"''. 

Finally, we note that comparison of the averages be- 
tween different declination bands are not meaningful; 
the above analysis is, therefore, insensitive to the ex- 
cess/deficit from the direction of the celestial poles. In 
other words, the 2-dimensional anisotropy can be thought 
of as a series of 1-dimensional curves in consecutive strips 
of declination. Before going through the analysis filter, 
each curve is described by a constant offset correspond- 
ing to the average flux, and a sum of harmonics whose 
average is zero. The filter removes the constant offset, 
keeping all other terms intact. Thus, the analysis pre- 
sented here is insensitive to any anisotropy along Earth's 
rotation axis. 



IV. THE ROBUSTNESS OF THE OBSERVED 
ANISOTROPY 

The result of the analysis is insensitive to the exact 
choice of the track length and zenith angle cuts. As an il- 
lustration of this inscnsitivity, track-type plots were made 
for 60 combinations of track length and zenith angle cuts; 
the former was varied between and 20 m, and the latter 
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TABLE I: Summary of one dimensional anisotropy measurements from deep underground muon telescopes. Col. 1: SK-l, KAM, 
and MAC refer to the SK-I, Kamiokande [T^], and MACRO [iJl experiments. Col. 2: type of plot. TRACK = track-type plot, 
ZENITH = zenith- type. CORR. = plot corrected for spurious sidereal anisotropy of atmospheric origin, UNC. = plot uncorrected 
for this. Col. 3: nominal value of detector projected area, in m^. Col. 4: nominal value of the overburden, in m.w.e.. Col. 
5: total detector live time, in days. Col. 6: total number (millions) of events. Cols. 7-10: Best fit first and second harmonic 
amplitude and phase. Errors are statistical except entries with two errors, where the first error is statistical and the second 
is the systematic error introduced in subtracting the atmospheric effect. Col. 11: P^r degree of freedom of fit of Egn. 1181 
to the data (x2 does not apply to data corrected for the atmospheric effect). The second harmonic amplitude and phase are 
the same for the corrected and uncorrected result because the spurious atmospheric anisotropy is assumed to vary as a first 
harmonic function. Kamiokande and MACRO report only a first harmonic fit to their data. 
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FIG. 4; Sky map of the anisotropy in equatorial coordinates. The sky is divided into 10° x 10° cells (Gouraud smoothing applied 
only for visual purposes). Declinations less than —53.58° (white region) always lie below the horizon and are thus invisible to 
the detector. In (a), each cell shows the fractional variation from the isotropic flux, while in (b) it shows the standard deviation 
of this variation. The solid red and blue curves show the excess and deficit cones obtained using a clustering algorithm applied 
to the data. The dashed curves in (a,b) show excess and deficit cones from the NFJ model Q, which is described in Section [V] 
(c) and (d) show the maps in (a) and (b) transformed to the galactic coordinates. The solid red and blue curves are the same 
cones as described above. The dashed red horizontal line indicates the direction of the Orion arm. The white patch indicates 
the below-horizon region. 
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(0.140 ±0.026)% 


TAIL-IN 


NFJ MODEL 


(90°, -24°) 


68" 




DEFICIT 


VIRGO 


OBSERVED, CORR. 


(205°, 5°) 


54° 


(-0.094 ± 0.014)% 


OBSERVED, UNC. 


(205°, 5°) 


54" 


(-0.099 ± 0.014)% 


GALACTIC 


NFJ MODEL 
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TABLE II: Cone parameters of excess and deficit regions. The NFJ model is described in Section [V] Column 3 describes 
whether the cone is from observation or model, and whether or not the atmospheric correction has been applied. Columns 4 
and 5 show the center and half opening angle of the cones. Column 6 shows the deviation from the isotropic event rate in the 
contained regions (the error is the statistical error based on the number of events contained in the cone). Rows labeled "CORR." 
and "uNC." refer, respectively, to cones obtained from the anisotropy map corrected and uncorrected for atmospheric effects. 



between 30° above horizontal to 90° below horizontal (i.e. 
no zenith angle cut). The harmonic function Egn.fTSlwas 
fit to each plot, and the RMS spread for each of the four 
parameters was found to be within 50% of the statistical 
error. 



As a test of signal robustness, the data were divided 
into five exact-year periods spanning June 1^* to May 31** 
of every year from 1996 to 2000, and a measurement of 
anisotropy from the track-type plot was made on each 
set. The best fit first harmonic amplitude and phase are 
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FIG. 5: The 2-parameter 68% confidence level regions {^y^ ~ 
2.3) of the amplitude and phase of the first harmonic function 
fit to yearly track-type anisotropy plots. The radial distance 
from the origin is the first harmonic amplitude, while the 
counterclockwise angle from RA = 0° is the right ascension 
at maximum. The regions are labeled by the corresponding 
year. The label "Combined" indicates the contour from the 
5-year combined data set. 

shown in Fig. [5l together with their 2-parameter 68% 
confidence level regions. Good overlap is seen. The fact 
that the phase is consistent from year to year is strong 
evidence that the observed anisotropy is due to a real 
physical effect. 



V. DISCUSSION 

The near-isotropy of cosmic rays with energy per nu- 
cleus below the "knee" in the spectrum is usually de- 
scribed as follows. At these energies, strong evidence 
exists to indicate that primary cosmic ray nuclei mostly 
originate as interstellar matter in the Milky Way Galaxy. 
They are accelerated by blast waves from supernova rem- 
nants, and are effectively trapped in the Milky Way 
by the Galactic magnetic field with a strength of order 
several micro-Gauss. The gyro-radius of a 10^^ eV ^ 
10^** eV proton propagating in a uniform micro-Gauss 
level magnetic field is in the range ~ 10 AU and ^ 0.1 pc, 
much smaller than the thickness of the Galactic disk of 
200 ~ 300 pc. The motion of cosmic ray nuclei is spiral- 
like in regions of the Galaxy where the magnetic field is 
smooth. Interspersed in these regions are areas where the 
magnetic field is irregular, in which the trajectories are 
complex and the direction before and after entrance in 
these areas is nearly random. Over large distance scales, 



the irregular regions can be thought of as scatterers of 
cosmic rays. As a result, the average motion of cos- 
mic rays in the Galaxy is expected to be highly random, 
which is consistent with the observed near-isotropy of the 
flux. 

The anisotropy at these energies is presumably due 
to a number of mechanisms. At the largest scales, the 
distribution of cosmic ray sources, the large-scale config- 
uration of the Galactic magnetic field, the distribution of 
magnetic field irregularities, and the location of the solar 
system in the Galaxy are all expected to be contributing 
factors. At smaller scales, the magnetic field configu- 
ration in the neighborhood of the solar system and the 
distribution of the nearest cosmic ray sources should con- 
tribute. At the smallest scales, the solar magnetic field 
can be ruled out at these energies, although there are sug- 
gestions that the heliosphere may play a role at energies 
around ~ 1 TeV (e.g. Q). Of course, the Compton- 
Getting effect is expected to produce a dipole anisotropy 
on top of all of the above. The relative importance of 
each of these effects is not known. 

As stated in section IIIH Earth-based cosmic ray 
anisotropy measurements at these energies require the 
application of a filter to the data in order to remove the 
large uncertainty in zenith angle dependence of the cos- 
mic ray flux; this filter also removes the axisymmetric 
component of the anisotropy along Earth's rotation axis. 
The part of the anisotropy that comes through this fil- 
ter is robustly established by several experiments. The 
earliest map of the large scale anisotropy, referred to as 
the NFJ models was made by Nagashima, Fujimoto, and 
Jacklyn by combining data from several different experi- 
ments in the northern and southern hemispheres j^] . The 
excess and deficit cones were obtained by interpolating 
between one dimensional anisotropy measurements made 
in several different declination strips. Because the data 
used were from very different detector types (shallow un- 
derground muon telescopes vs. surface air shower array) 
with correspondingly large spread in energy sensitivity 
and very different systematic uncertainties, the result was 
qualitative in nature. More recently, large detectors with 
correspondingly large statistics and good pointing accu- 
racy have come on line, and each one is able to make 
a map of the large scale anisotropy. We, as well as the 
Tibet air shower experiment (27j . have published such 
maps, and both agree well with each other, as well as 
with that of [i]. 

A unique interpretation of the observed anisotropy is 
not possible, but it is useful to categorize the interpre- 
tations into two classes: (1) the true anisotropy is dom- 
inated by the dipole term, and (2) the higher harmonics 
are not negligible. If scenario (1) were true, then the dis- 
tortion introduced by the analysis method projects the 
dipole onto the equatorial plane. Also, an excess cone and 
deficit cone should exist in 180° opposition to each other. 
The cones found in this analysis are both centered close 
to the equatorial plane; however, they are separated in 
right ascension by 130°, which is in apparent contradic- 
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tion to the dipole-dominant hypothesis. Quantitatively, 
a fit of the data to an equatorial dipole gives a 13% 
chance that the observed map is consistent with a pure 
equatorial dipole hypothesis. Thus the data do not rule 
out this scenario. 

Before discussing scenario (2), let us consider the im- 
plications of this scenario. Considering the complicated 
nature of the origins and propagation of cosmic rays, it is 
reasonable to assume that several different mechanisms 
contribute to the overall observed anisotropy; each mech- 
anism contributes a dipole term to the map, and the ob- 
served dipole is a sum over all the dipoles projected onto 
the equatorial plane. One component of the dipole is due 
to the Compton-Getting effect. Since the nature of the 
other mechanisms is unknown, it is not possible to extract 
the Compton-Getting term. If, however, we make the ex- 
treme interpretation that the Compton-Getting effect is 
the dominant term, then one can extract the equatorial 
projection of the relative velocity between the cosmic ray 
rest frame and the solar system. The fit described 
above gives D = {7.lt]:l) ^ lO"** and uq = 32° ± 12° 
with xVd-O.f. = 577/538 (see Eqn. This corre- 

sponds to a relative velocity of 49l|^|!) km/s in the di- 
rection 32° ±12° right ascension. The speed is signifi- 
cantly smaller than the orbital speed of the solar system 
around the Galaxy (« 200 km/s), while it is comparable 
to the relative speed of neighboring stars around the sun. 
Unless there is an accidental cancellation of large dipole 
terms, the observed speed should be about the magni- 
tude of the actual Compton-Getting speed. Thus one 
can deduce that, very likely, cosmic rays in the neighbor- 
hood of the solar system move around the galaxy with a 
motion similar to stars. In other words, the cosmic ray 
rest frame is dragged along with stars. 

Let us now consider scenario (2), i.e. higher harmonic 
terms are not negligible. For this scenario, we focus on 
the particular form where two independent cones - one 
with an excess flux and the other with a deficit - produce 
the observed anisotropy. The right ascension of the cone 
center is not affected by the distortion, whereas the decli- 
nation may or may not be significantly distorted. Specif- 
ically, if the true declination of the excess cone center 
is similar to that of the deficit cone, then the observed 
declination value should be equal to the true value (to 
within statistical uncertainty). On the other hand, if 
there is a mismatch in the declination values, then the 
true values are farther apart than observed - i.e. the 
filter causes the reconstructed cone declination values to 
get pushed towards each other. According to Table [TTl 
the observed declination of the excess cone center is —5°, 
while that of the deficit cone is -1-5°; these two values 
are close to each other, indicating that the true values 
couldn't have been too far from the equator. To put 
this statement on a quantitative footing, the data were 
compared with anisotropy maps formed with different 
assumptions regarding the true parameter values of the 
excess and deficit cones. Each cone is defined by four 
parameters: the position of the cone center (two param- 
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FIG. 6: A map of the marginalized in the parameter space 
defined by the declination of the excess and deficit cones. At 
each point in the map, was minimized with respect to the 
four remaining parameters. 



eters), the opening angle, and the amplitude (assumed 
to be constant within the cone). The right ascension of 
each cone was fixed to the observed value since it is not 
distorted by the analysis filter. There are, thus, six free 
parameters that describe the two cones. Each set of pa- 
rameter values gives rise to a value of when compared 
against the data. The absolute minimum x^/d.o.f. of 
544/534 occurs near the set of values given in Table [Til 
A marginalized x^ map is shown in Fig. [SI The figure 
shows as a function of declination of the excess and 
deficit cone center. At each point in the map, x^ was 
minimized with respect to the remaining four parameters. 
The contour shows the confidence level with which any 
pair of declination values is allowed. At 90% confidence 
level, the declination of the excess cone is completely un- 
constrained. A more stringent constraint is attainable 
with the deficit cone declination, but at 99% level, it also 
becomes almost totally unconstrained. 

Let us explore the implications of the scenario where 
the observed cones are close to the true ones. A natural 
coordinate system for interpreting cosmic ray anisotropy 
in these energies is galactic. Figs. [J (c) and (d) show the 
anisotropy map in galactic coordinates, (c) showing the 
fractional variation from isotropy, while (d) shows the 
standard deviation of this variation. One feature that 
clearly stands out is the large deficit seen in the Galactic 
northern hemisphere. This may be related to the fact 
that the solar system is displaced to the north of the 
Galactic equator by about 15 pc, a significant fraction of 
the half-width of the Galactic disk of about 100^^200 pc. 
Since the cosmic ray density is almost certainly great- 



11 



est at the equator and tapering off with vertical distance 
away from it, the density is less to the north than to the 
south as viewed from the solar system. This will tend to 
produce a deficit flux in the Galactic north. Another ob- 
servation in this regard is the fact that the deficit is less 
pronounced in the direction of the Orion arm. The solar 
system is currently located at the edge of this arm and 
moving towards it. It is a standard view that the cosmic 
ray density is elevated in the spiral arms compared to 
the gap regions. Thus a density gradient is likely to pro- 
duce an excess flux from the Orion arm as viewed from 
the solar system; this excess flux may be canceling out 
the deficit flux from the Galactic north. We note that 
these observations essentially the same as those made 
in [2^. A final unexplained feature is the excess cone. 
We are unaware of any Galactic features that may cause 
this. The proponents of the NFJ model Q point out 
that it is more-or-less aligned with the tail direction of 
the heliosphere; however, no plausible physical mecha- 
nism exists that could explain the size of the observed 
anisotropy poj . 



VI. CONCLUSION 

An anisotropy map of cosmic rays of nominal en- 
ergy of 10 TeV was made from 1662 days of obser- 
vation. The right ascension projection of this map 
has a first harmonic amplitude and phase of (6.64 ± 
0.98 (stat.)±0.55 (syst.))xlO-* and (33.2°±8.2° (stat.)± 
5.1° (syst.))°, which are in good agreement with results 
from other experiments. The sky map indicates a region 
with (0.104 ± 0.020)% excess flux in the constellation of 
Taurus, while a region with (0.094±0.014)% deficit flux is 
observed in the constellation of Virgo. The excess region 
is centered at (ax, St) = (75° ± 7°, -5° ± 9°) with a half 
opening angle of 39° ± 7°, while the corresponding values 
for the deficit region are {av,Sv) (205° ± 7°, 5° ± 10°) 
and half opening angle = 54° ± 7° . These regions largely 
coincide with those of the NFJ model, and also with those 
observed by the Tibet collaboration. This agreement 
between experiments using very different measurement 
techniques spanning several decades of observation and 
covering primary cosmic ray energies of ~ 1 to ~ 100 TeV 
indicates the robustness of the observed anisotropy pat- 
tern. The pattern, therefore, is a real feature of the cos- 
mic ray flux in the neighborhood of the solar system at 
the current epoch. The simplest model for the observa- 
tion is a pure dipole pattern, which could be produced by 
the Compton-Getting effect, but could also be the lead- 
ing harmonic term from other more complicated mecha- 
nisms. Our observation is not described very well with 
a pure dipole pattern, although, at 13% confidence level, 
we cannot rule it out. The distorting effect of the filter 
applied to the data prevents us from making a unique 
physical interpretation of the observation. If, however, 
it is assumed that the distortion is not too great, the 
deficit region coincides with a large portion of the Galac- 



tic northern hemisphere. This may be related to the fact 
that the solar system is displaced by about 15 pc to the 
north of the Galactic plane. Also, the deficit appears to 
be canceled out in the direction of the Orion Arm, which 
may be an indication of enhanced levels of cosmic ray 
density there. The excess region does not seem to match 
any features that could provide a mechanism for its exis- 
tence, though it has been noted [1] that it points in the 
direction of the tail end of the heliosphere. 
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APPENDIX A: SUBTRACTING THE 
ANISOTROPY DUE TO ATMOSPHERIC 
EFFECTS 

The observed sidereal anisotropy is due to two effects: 
extra-terrestrial (i.e. Compton-Getting, Galactic, and 
heliospheric effects) and atmospheric (due to residual ef- 
fects of seasonal and solar diurnal variations in the atmo- 
spheric temperature) . Since we are interested only in the 
extra-terrestrial anisotropy, it is desirable to subtract the 
atmospheric contribution. We discuss in this Appendix 
technical details on this subtraction. The method is due 
to Farley and Storey [2], which is applied to the zenith- 
type plot. In the second section this result is general- 
ized to the two dimensional anisotropy map. Finally, the 
results of the first two sections are used to subtract the 
atmospheric contribution from the track-type one dimen- 
sional plot. 



1. Subtraction for the Zenith- Type Plot 

The zenith-type plot is equivalent to the relative varia- 
tion in the muon rate as a function of local sidereal time. 
In Farley and Storey's formulation, the rate variation is 
parameterized generally as follows: 
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SEASONAL MODULATION SOLAR VARIATION 

, V , ^ 

R{t) = 1 + [A + 2B cos27r(< - (?!)2)] cos 27r(A^t - 0i) 

+ C cos2^{(Af + l)i-03}, (Al) 
. ' 

TRUE SIDEREAL VARIATION 

where t is measured in years, N « 365.24 cycles/year is 



the solar diurnal frequency, and c/ii, i = 1,2,3 are phase 
angles. The parameters A, B, and C are the magnitude 
of the relative rate variation for different periodicities 
(discussed below). The solar diurnal variation is assumed 
to be seasonally modulated (first line, Eqn. lAip . The 
second line represents the true sidereal variation (i.e. of 
extra-terrestrial origin). A re-arrangement of the first 
line above gives the following: 



SOLAR 



R{t) = l + Acos2Tr{Nt-(f)i) (A2) 

SPURIOUS SIDEREAL TRUE SIDEREAL 



+ B' cos2tt{{N + l)t~ {(f>i +(p2)} + C cos2tt{{N +l)t- <j>3} 
+ B COS 27r {(iV - l)t - {(jji - (/)2)} 



PSEUDO-SIDEREAL 



The first time dependent term in Eqn. IA2I is the rela- 
tive rate variation with a periodicity of one solar day, 
the second and third terms are rate variations with a pe- 
riodicity of one sidereal day, and the final term is the 
rate variation with a periodicity of one pseudo-sidereal 
day, which is longer than the solar diurnal day by about 
0.27% (Figs. [7| (a)-(c)). Written in this form, it is seen 
that there are two sources of sidereal variation, the "spu- 
rious" one due to the atmosphere and the "true" one due 
to extraterrestrial effects. 

Each time dependent term in Eqn. IA2I can be repre- 
sented by phasors. In Cartesian coordinates, they are 
given as follows: 



A 
B' 
C 
B 



= {A cos( 



ii, A sin0i) 
{B' cos(0i -1-02), B' sin(<^i -1-02)) 
(C cos 03, C sin 03) 
{B cos(0i - 02), B sin(0i - 02)) 

The phasor A is non-zero due to residual effects of the 
solar diurnal and seasonal variation in the atmospheric 
temperature, while D = C + B' is non-zero primarily 
due to extra-terrestrial effects (represented by C), al- 
though, in general, a non-zero contribution is also made 
by atmospheric effects (represented by B'). No real ef- 
fect is directly responsible for a non-zero value of B, but, 
as described above, interplay between seasonal and so- 
lar diurnal variation in the atmospheric temperature can 
indirectly give rise to a non-zero magnitude. 

In terms of phasors, one sees that the process of mea- 
suring the true sidereal variation involves measuring the 
phasor B' and subtracting it from D. The phasor D is 
obtained from Fig. [7] (a), while B' is obtained from B 
and A, which are, in turn, obtained from Figs. [7] (b) and 
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FIG. 7: (a) Relative muon rate as a function of local sidereal 
time (i.e. zenith-type map), in hours right ascension, (b) Rel- 
ative muon rate as a function of local solar hour, (c) Relative 
muon rate as a function of hours, pseudo-sidereal time. The 
curve in each frame is the best fit of the first two harmonic 
functions to the data. 



(c) . Specifically, B' is obtained by reflecting B about the 
axis defined by A (see Fig. |S] (a)). 

B' and C can be obtained approximately by using the 
most likely value of A, B, and D. These are given as 
i3'(CALC) and (7(CALC) in Table Hill A statistically rig- 
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FIG. 8: Phasor diagrams showing the result of subtracting the atmospheric and Compton- Getting effects. The length of an 
arrow represents the amplitude of the first harmonic component, while the angle measured counter clockwise from <j) — 0° 
is the phase at maximum, (a) Atmospheric effect, (b) Compton- Getting effect assuming cosmic ray rest frame moving with 
the local standard of rest, and (c) same as (b), but moving with the local interstellar medium. In each plot, the vector D 
indicates the uncorrected amplitude and phase^ while the vector —B', — Vlsr, and — Vism are corrections for the atmospheric 
and Compton-Getting effects. The vectors A, B, and C are defined in the text, (b) and (c) also show the result of subtracting 
both effects. 



PHASOR 


AMP (xlO ") 


PHASE (hr) 


A 


6.63 ±0.98 


18.0 ±0.6 


B 


1.01 ±0.98 


13.1 ±3.7 


D 


5.65 ±0.98 


2.3 ±0.7 


B' (CALC) 


1.01 


22.9 


C (CALC) 


5.09 


2.9 


B' (stat) 


0.98 ±0.71 


23.2 ±5.2 


C (stat) 


5.31 ± 1.19 


2.7 ±0.9 



TABLE III: Summary of phasor parameters. A, B, and D are 
from data (N.B. the amplitude and phase of D are slightly 
different from those shown in Table |I] because of difference in 
binning); i3'(CALC) and (7(CALC) are calculated directly from 
the first three, while B'(stat) and (7(stat) are obtained by 
the statistical subtraction technique described in the text. 



orous determination of C and associated uncertainties 
require the use of a statistical subtraction technique in 
which the phasors A, B, and D are generated randomly 
according to their respective probability from the first 
harmonic fits in Fig.[71 For each generated triplet of pha- 
sors, there corresponds a unique value of C, and an en- 



C 



and ( 



re- 



semble of generated C gives a distribution of 

The mean and RMS of a Gaussian fit to each are taken as 
the most likely value and uncertainty of the true sidereal 
anisotropy; these are given as (7(stat) in Table Hill 



dimensional map, the true (i.e. corrected) map A^TRUE IS 
given by: 



Mt 



Mo 



Ma 



(A3) 



where Mobs is the observed map, and Matm is the map 
of anisotropy of atmospheric origin. Matm is calculated 
by assuming that: (1) the incident cosmic ray fiux is 
isotropic (the observed anisotropy introduces only a sec- 
ond order correction, so it can be ignored); (2) the atmo- 
spheric effect causes the overall cosmic ray rate to vary 
with amplitude and phase given by the phasor _B', i.e. 
amplitude = 0.98 x 10~* and phase = 348° right ascen- 
sion. 

The map Matm is obtained by convoluting the relative 
rate variation R{ts) = 1 + \B'\ cos(a;s Ts — (pB') with the 
isotropic event rate I{S,h), where Ts is sidereal time, ujs 
is the sidereal angular frequency, S is the declination, 
h — a ^ Ts is the hour angle, and a is the right ascension. 
Note that / has units of day~^ m~^ sr~^; it is related to 
Fig. [T] by a coordinate transformation. The convolution 
is as follows: 



Matm(i5, a) 



J dTsI{6,a 



UsTs) R{Ts) 



(A4) 



2. Subtraction for the Two Dimensional Map 

We describe in this section the method used to sub- 
tract the anisotropy of atmospheric origin from the two 
dimensional anisotropy map. Letting M denote a two 



The map Matm (from which the S dependence is factored 
out) is shown in Fig. [3] (a). The excess and deficit cone 
parameters for Mtrue are shown in Table IIVI The direc- 
tion and cone size of the deficit region are unchanged by 
this correction, whereas those of the excess region change 
noticeably. 



14 



REGION TYPE 


NAME 


CONE SOURCE 


(Q,d) 


SIZE 


DEVIATION 


X 


EXCESS 


TAURUS 


NO CORRECTION 


(65°, 5°) 


27° 


0.140% 


5.26 0- 


ATM 


(75°, "5") 


39° 


0.104% 


5.31(7 


CG, LSR 


(85°, -35°) 


65° 


0.0964% 


6.46(7 


CG, ISM 


(55 ,5 ) 


6U 


0.088470 


0.90 (7 


ATM + CG, LSR 


(80 , — 3o ) 


60 


o.oyyiTo 


0.00 (7 


ATM + CG, ISM 


(00 , j 


6U 


O.Oe 1 I/O 


0.9U (7 


TAIL-IN 


NFJ MODEL 


(90 ,"24 ) 


68° 






DEFICIT 


VIRGO 


NO CORRECTION 


(205°, 5°) 


54° 


-0.0988% 


-7.27(7 


ATM 


(205°, 5°) 


54° 


-0.0940% 


-6.91(7 


CG, LSR 


(215 ,5 ) 


00 


— 0.iU( 70 


— /.90 (7 


CG, ISM 


(215°, 5°) 


55° 


-0.115% 


-8.57a 


ATM + CG, LSR 


(215°, 5°) 


55° 


-0.103% 


-7.66(7 


ATM + CG, ISM 


(215°, 5°) 


55° 


-0.111% 


-8.28(7 


GALACTIC 


NFJ MODEL 


(180°, 20") 


57° 







TABLE IV: Cone parameters of excess and deficit regions. Tlie column "Cone Source" refers to the data set or model from 
which the cones are derived. "No Correction" refers to the anisotropy without any subtraction, while the other rows refer 
to that after subtracting various eflects. "ATM" is the atmospheric effect, "CG, LSR (ISM)" is the equatorially projected 
Compton-Getting effect with the cosmic ray rest frame assumed to be the same as the local standard of rest (local interstellar 
matter). "ATM -|- CG, LSR (ISM)" is that from which both effects are subtracted. Columns 4 and 5 show the center and half 
opening angle of the cones. Column 6 shows the deviation from the isotropic event rate in the contained regions, and column 
7 (x) shows the statistical significance of the deviation. 



3. Subtraction for the Track-Type Map 

The track-type one dimensional map corrected for at- 
mospheric effects can be obtained by simply projecting 
Mtrue((5, a) onto the right ascension axis. However, this 
does not provide an estimate of the uncertainty intro- 
duced by the atmospheric subtraction. In order to ob- 
tain this, we first generate an ensemble of phasors D, A, 
and B, as described in Sec. lA II For each triplet, a two 
dimensional map Mtrue is made, as described in Sec IA 21 
This map is then projected onto right ascension axis to 
obtain the one dimensional map. The first two harmonic 
functions (Eqn. [T5|) are then fit to each map thus ob- 
tained. We thus obtain an ensemble of fit values Ai and 
4>i (the second harmonic is unchanged in the ensemble 
because the atmospheric effect was assumed to have only 
first harmonic variation). The result of this procedure is 
given in Table |T] in the row labeled "track/corr.". The 
first error is statistical (from fitting harmonic functions 
to the data), and the second error is from the disper- 
sion in the fit values obtained from the ensemble method 
described above. 



APPENDIX B: CORRELATION BETWEEN 
MUON RATE AND ATMOSPHERIC 
TEMPERATURE 

The effect of the atmosphere on the cosmic ray detec- 
tion rate in underground muon detectors is, in general, 
correlated with the pressure at the detector altitude and 
with the atmospheric temperature profile above the de- 
tector. The pressure dependence becomes unimportant 
compared to temperature dependence for muon threshold 



energy greater than about 100 GeV [sot . 

In this limit, the relative change in the muon rate with 
atmospheric temperature is given by the following ex- 
pression: 

61 r° — 

— « / a{x, Eo,xo) ST{x) dx (Bl) 
Jq 

The quantity a is the partial temperature coefficient, 
6T{x) is the deviation of the temperature from the mean 
at atmospheric depth x, Eq is the threshold muon en- 
ergy, and xq is the atmospheric depth at the detector al- 
titude. The mechanism for the temperature dependence 
of the rate is as follows. As the temperature rises, the 
atmospheric density decreases, and the probability that 
a meson in a cosmic ray shower is destroyed by interac- 
tion with air nuclei decreases. The increased meson mean 
free path implies that mesons have increased chance to 
decay and produce muons. The cosmic ray muon rate, 
therefore, is positively correlated with atmospheric tem- 
perature. The partial temperature coefficient can be cal- 
culated numerically using inputs such as a model atmo- 
sphere, primary cosmic ray flux, particle production cross 
section, particle decay constants, etc.l30|, while 6T{x) 
can be obtained at discrete atmospheric levels from me- 
teorological measurements. For SK-I, the change in rate 
due to this effect should be « ±1%, which is more than 
an order of magnitude larger than the magnitude of the 
sidereal anisotropy. 

Figure [TO] (a) shows the relative variation in the muon 
rate for each month and year of SK-I. The solid curve 
shows the predicted variation based on Egn. lBll 31]. The 
temperature measurements were obtained from the Wa- 
jima Observatory (37.38° N, 136.90° E, 116 km from the 
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FIG. 9: Anisotropy introduced by: (a) the atmospheric effect; 
(b) the Compton-Getting effect assuming that the bulfc cosmic 
ray motion is the same as the local standard of rest; (c) same 
as (b), but the motion is assumed to be the same as that of 
the neutral interstellar matter. The contour values indicate 
fractional deviation from isotropy. The white region below 
declination of —53.58° is always below the horizon. Note that 
the filter applied to the data projects the anisotropy shown 
in (b) and (c) onto the equatorial plane. 



spheric depths, so inaccurate temperature measurements 
at altitude above 40 mb (i.e. pressure < 40 mb) introduce 
significant inaccuracies in the predicted rate. 

A better standard for checking the SK-I rate variation 
is by means of a simultaneous independent measurement 
from a nearby underground muon detector. The Mat- 
sushiro detector (36.53 N, 138.01 E, 79 km from SK-I) is 
perfectly suited for this requirement [33| . With an over- 
burden of 220 m.w.e., its muon energy threshold is about 
100 GeV. This lower energy threshold implies that the 
rate variation at Matsushiro is similar to that of SK-I, 
but has smaller amplitude (see Fig. [TO] (b)). 

One measure of this difference is the month to month 
variation in AJ// at the two sites. According to calcula- 
tion 31], the magnitude of this change at SK-I should be 
2.03 times larger than that at Matsushiro. A correlation 
plot of the changes at the two sites is shown in Fig. [11] 
The regression coefficient /? = 2.03 ±0.05 is in agreement 
with the predicted value. 

When the data are binned in solar diurnal hours, the 
1% level monthly variations almost cancel out, leav- 
ing a residual variation at the level of several parts 
per ten thousand. This variation, when modulated sea- 
sonally, produces side band components with frequency 
365.24 ± 1 cycles per year. The frequency of 366.24 cy- 
cles per year is the inverse of one sidereal day, and the 
existence of this component implies that the observed 
sidereal variation in the cosmic ray rate is partly due to 
atmospheric temperature variations. This contribution 
to the observed sidereal anisotropy can be estimated us- 
ing the method of Farley and Storey . A detailed dis- 
cussion of the atmospheric contribution to the sidereal 
anisotropy is given in Appendix [^ 



APPENDIX C: SUBTRACTING THE 
COMPTON-GETTING ANISOTROPY 

The Compton-Getting effect refers to the enhancement 
of the cosmic ray flux in the observer's direction of motion 
relative to the reference frame in which the bulk motion 
of the cosmic ray plasma is at rest. If the observer's ve- 
locity relative to the cosmic ray bulk motion is v and the 
direction of observation is in the direction of the unit vec- 
tor M, the relative enhancement in the intensity is given 
by: 



— = 2 + 7 - cosx, 
1 c 



(CI) 



SK-I detector) [S^l. Radio sonde was used to measure 
the temperature of 25 layers of the atmosphere between 
1000 mb and 5 mb. These measurements were made twice 
a day. The agreement between the data and prediction is 
good, though not perfect. The disagreement is due to in- 
accuracies in the temperature measurements at altitudes 
above 40 mb. At SK-I energies, the partial temperature 
coefficient increases all the way to very shallow atmo- 



where v = \v\, c is the speed of light in vacuum, and 
cosx = V ■ u/v is the cosine of the opening angle be- 
tween the observer's motion and the direction of obser- 
vation. The cosmic ray rest frame is not known, al- 
though the observed smallness of the anisotropy implies 
that its motion relative to the sun must be small {v/c 
must be on the order of 10~^, or w < 30 km/s, barring 
an accidental large cancellation of the Compton-Getting 
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FIG. 10: (a) Variation in the muon rate relative to the mean for each month of SK-I. (b) The variation seen at Matsushiro 
during the same period. The solid curve in each frame is the predicted rate variation based on equation IBll 



anisotropy by that due to other effects). Two assump- 
tions are often invoked in the literature: the cosmic ray 
rest frame is at rest with respect to (1) the local stan- 
dard of rest, or (2) the local interstellar medium. The 
first motion has speed Ulsr ~ 20 km/s in the direction 
(aLSR, <^lsr) ~ (270°, 29.2°), while the second motion has 
speed VisM ~ 22 km/s in the direction (aisM, <5ism) ~ 
(252°, —17°). The values were chosen to be consistent 
with Q ; see references therein for citations for these val- 
ues. The dipole anisotropy due to the Compton-Getting 
effect is shown in Figs. |9] (b) and (c). As mentioned in 



Section IIIIl the filter applied to the data projects the 
dipole onto the equatorial plane, so the center of the ob- 
served dipole has no declination component, and the ef- 
fective velocity is the equatorial projection of the velocity, 
which is i;lsr cos i5lsr = 17 km/s for motion with respect 
to the local standard of rest, and ?;ismCOsJism = 21 km/s 
for motion with respect to the interstellar medium. The 
cone parameters before and after subtracting the equato- 
rial projection of the Compton-Getting anisotropy from 
the observed anisotropy are summarized in Table [TVl 
We also examined the effect of the Compton-Getting 
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FIG. 11: The correlation between (A///)„,„„i, SK-1 t;s. 
Matsushiro. AI/I is the relative rate variation of a given 
month; the notation (A7/7)„,„_i indicates the difference of 
this quantity between month n and n — 1. 7 = 0.855 is the 
linear correlation coefficient between the two data (60 degrees 
of freedom), while (3 = 2.03 ± 0.05 is the best fit slope. 



anisotropy on the one dimensional anisotropy. Since 
the track- and zenith-type plots are very similar, we 
focus here on just the zenith-type plot. Starting with 
the two dimensional Compton-Getting anisotropy map 
(Figs. [9] (b) and (c)) and folding in the effect of the 
overburden (Fig.[T]), the one dimensional anisotropy due 
to the Compton-Getting effect alone is described well 
by a first harmonic function with amplitude and phase 
(1.66 X 10"'*, 274°) for cosmic ray bulk motion with the 
local standard of rest, and (2.00x 10""*, 256°) for bulk mo- 
tion with the interstellar matter. These two anisotropics 
are indicated by the phasors I^sr and V^sm in Fig. [8] (b) 
and (c). The amplitude and phase of the observed 
anisotropy without any subtraction are indicated in the 
figure as D. The contribution of the Compton-Getting 
anisotropy is removed by subtracting ^LSRasM from D: the 
result is also shown in the figure. The anisotropy is en- 
hanced, and the direction of anisotropy rotates toward 
90°. The effect of atmospheric subtraction is also shown. 
The amplitude and phase of the anisotropy before and af- 
ter subtraction are summarized in Table IVl Also shown 



in the table are the anisotropics after subtracting both 
effects. 



APPENDIX D: THE CLUSTERING ALGORITHM 



The clustering algorithm is applied to a histogram Ni_j , 
where Nij is the exposure-corrected number of events in 
bin 



SUBTRACTION 


AMPLITUDE 


PHASE (dec) 


NONE 


5.7 X 10"* 


35° 


CO, LSR 


6.7 X 10"* 


48° 


CG, ISM 


7.3 X 10"* 


45° 


ATMOSPHERIC 


5.3 X 10-* 


40° 


CG, LSR -I- ATMOSPHERIC 


6.5 X 10-* 


52° 


CG, ISM -1- ATMOSPHERIC 


7.0 X 10-* 


50° 



TABLE V: Amplitude and phase of the zenith-type one di- 
mensional anisotropy before and after subtractions. CG, LSR 
refers to the Compton-Getting effect assuming cosmic ray rest 
frame moving with the local standard of rest, while CG, ISM 
refers to the case where the cosmic ray is assumed to move 
with the local interstellar matter. 



bin index. For each (i,j), the quantity x is calculated 
over a variable sized cone centered on (i, j), where: 



(Dl) 



'hs 



is the observed number of events in the cone, and 
is the expected number in the absence of anisotropy. 
The cone size (half opening angle) that extremizes x is 
sought; the excess/deficit is assumed significant if |x| > 4. 
If the center of one cone falls within another cone, the 
cone with small |x| is rejected. 

The statistical error of the reconstructed cone param- 
eters was estimated using an ensemble experiment tech- 
nique in which an input sky map with anisotropy cones 
with parameters given in Table [IT] was used to generate 
a large number of output maps with random statistical 
fluctuations. The clustering algorithm was applied to 
each generated map, and the distribution of the recon- 
structed cone parameters was examined. The RMS of 
these distributions were taken as the statistical error. 
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